{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4dab0844",
   "metadata": {},
   "source": [
    "# Tutorial 04: first eigenvalue of the Laplacian with Dirichlet BCs\n",
    "\n",
    "In this tutorial we evaluate the first eigenvalue of the Laplacian with homogeneous Dirichlet boundary conditions. Let $\\Omega$ be the unit ball in 2D: given a constant $\\alpha \\in \\mathbb{R}^+$, the goal is the find the smallest eigenvalue $\\eta \\in \\mathbb{R}^+$ such that\n",
    "$$\\begin{cases}\n",
    "-\\alpha \\Delta u = \\eta u, & \\text{in } \\Omega,\\\\\n",
    " u   = 0, & \\text{on } \\partial\\Omega,\n",
    "\\end{cases}$$\n",
    "\n",
    "The weak formulation of the eigenvalue problem is\n",
    "\\begin{align*}\n",
    "&\\text{find } (\\eta, u) \\in \\mathbb{R} \\times V \\text{ s.t. }&\\\\\n",
    "&\\alpha \\int_\\Omega \\nabla u \\cdot \\nabla v = \\eta \\int_\\Omega \\; u \\; v, & \\forall v \\in V,\\\\\n",
    "\\end{align*}\n",
    "where\n",
    "$$\n",
    "V = H^1_0(\\Omega).\n",
    "$$\n",
    "\n",
    "In the following, we will adopt the notation $\\eta = \\eta^{(\\alpha)}$ when interested in comparing the first eigenvalue for different values of the parameter $\\alpha$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85009d8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import typing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8f1a4aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "import dolfinx.fem\n",
    "import dolfinx.fem.petsc\n",
    "import dolfinx.io\n",
    "import gmsh\n",
    "import mpi4py.MPI\n",
    "import numpy as np\n",
    "import petsc4py.PETSc\n",
    "import scipy.special\n",
    "import slepc4py.SLEPc\n",
    "import ufl\n",
    "import viskex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad2e0524",
   "metadata": {},
   "outputs": [],
   "source": [
    "import multiphenicsx.fem\n",
    "import multiphenicsx.fem.petsc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1d83944",
   "metadata": {},
   "source": [
    "### Geometrical parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84e72abc",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh_size = 0.05"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "092d4b9e",
   "metadata": {},
   "source": [
    "### Mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cd7acb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "gmsh.initialize()\n",
    "gmsh.model.add(\"mesh\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d7793b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "p0 = gmsh.model.geo.addPoint(0.0, 0.0, 0.0, mesh_size)\n",
    "p1 = gmsh.model.geo.addPoint(0.0, 1.0, 0.0, mesh_size)\n",
    "p2 = gmsh.model.geo.addPoint(0.0, -1.0, 0.0, mesh_size)\n",
    "c0 = gmsh.model.geo.addCircleArc(p1, p0, p2)\n",
    "c1 = gmsh.model.geo.addCircleArc(p2, p0, p1)\n",
    "boundary = gmsh.model.geo.addCurveLoop([c0, c1])\n",
    "domain = gmsh.model.geo.addPlaneSurface([boundary])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1c51a18",
   "metadata": {},
   "outputs": [],
   "source": [
    "gmsh.model.geo.synchronize()\n",
    "gmsh.model.addPhysicalGroup(1, [c0, c1], 1)\n",
    "gmsh.model.addPhysicalGroup(2, [boundary], 0)\n",
    "gmsh.model.mesh.generate(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e19c55b",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh, subdomains, boundaries = dolfinx.io.gmshio.model_to_mesh(\n",
    "    gmsh.model, comm=mpi4py.MPI.COMM_WORLD, rank=0, gdim=2)\n",
    "gmsh.finalize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d81ab434-9b59-4279-80e5-aeea6807a3f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create connectivities required by the rest of the code\n",
    "mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb1744ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "facets_partial_Omega = boundaries.indices[boundaries.values == 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40817c9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "viskex.dolfinx.plot_mesh(mesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37ed37ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "viskex.dolfinx.plot_mesh_tags(mesh, boundaries, \"boundaries\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f98dd8f",
   "metadata": {},
   "source": [
    "### Approach #1\n",
    "\n",
    "We define a $\\mathbb{P}^2$ finite element space. In the simplest approach, we then use `dolfinx` to assemble the left-hand side matrix $A = A^{(\\alpha)}$ representing the discretization of the Laplace operator, and the right-hand side matrix $B$ representing the discretization of the $L^2(\\Omega)$ inner product. We then apply Dirichlet boundary conditions through `dolfinx`, and compare the obtained first eigenvalue with the analytical solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34f88662",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a function space\n",
    "V = dolfinx.fem.functionspace(mesh, (\"Lagrange\", 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e54e15b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define trial and test functions\n",
    "u = ufl.TrialFunction(V)\n",
    "v = ufl.TestFunction(V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4178408b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Auxiliary function for the computation of the first eigenvalue using dolfinx\n",
    "def get_smallest_eigenvalue_dolfinx(  # type: ignore[no-any-unimported]\n",
    "    alpha: petsc4py.PETSc.RealType, diagonal_A: petsc4py.PETSc.RealType = 1.0\n",
    ") -> tuple[\n",
    "    petsc4py.PETSc.RealType, dolfinx.fem.Function\n",
    "]:\n",
    "    \"\"\"Get the smallest eigenvalue, and the corresponding eigenfunction, using dolfinx.\"\"\"\n",
    "    # Define problem forms\n",
    "    alpha_constant = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(alpha))\n",
    "    a = alpha_constant * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx\n",
    "    b = ufl.inner(u, v) * ufl.dx\n",
    "    # Define boundary conditions\n",
    "    zero = petsc4py.PETSc.ScalarType(0)\n",
    "    bdofs_V = dolfinx.fem.locate_dofs_topological(V, mesh.topology.dim - 1, facets_partial_Omega)\n",
    "    bc = dolfinx.fem.dirichletbc(zero, bdofs_V, V)\n",
    "    # Assemble lhs and rhs matrices\n",
    "    A = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(a), bcs=[bc], diagonal=diagonal_A)\n",
    "    A.assemble()\n",
    "    B = dolfinx.fem.petsc.assemble_matrix(dolfinx.fem.form(b), bcs=[bc])\n",
    "    B.assemble()\n",
    "    # Solve\n",
    "    eps = slepc4py.SLEPc.EPS().create(mesh.comm)\n",
    "    eps.setOperators(A, B)\n",
    "    eps.setProblemType(slepc4py.SLEPc.EPS.ProblemType.GHEP)\n",
    "    eps.setDimensions(1, petsc4py.PETSc.DECIDE, petsc4py.PETSc.DECIDE)\n",
    "    eps.setWhichEigenpairs(slepc4py.SLEPc.EPS.Which.SMALLEST_REAL)\n",
    "    eps.getST().getKSP().setType(\"preonly\")\n",
    "    eps.getST().getKSP().getPC().setType(\"lu\")\n",
    "    eps.getST().getKSP().getPC().setFactorSolverType(\"mumps\")\n",
    "    eps.solve()\n",
    "    assert eps.getConverged() >= 1\n",
    "    # Extract first eigenvalue and eigenvector\n",
    "    vr = dolfinx.la.create_petsc_vector(V.dofmap.index_map, V.dofmap.index_map_bs)\n",
    "    vi = dolfinx.la.create_petsc_vector(V.dofmap.index_map, V.dofmap.index_map_bs)\n",
    "    eigv = eps.getEigenpair(0, vr, vi)\n",
    "    r, i = eigv.real, eigv.imag\n",
    "    assert abs(i) < 1.e-10\n",
    "    assert r > 0., \"r = \" + str(r) + \" is not positive\"\n",
    "    # Destroy EPS object\n",
    "    eps.destroy()\n",
    "    # Transform eigenvector into an eigenfunction that can be plotted\n",
    "    vr.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)\n",
    "    vr_fun = dolfinx.fem.Function(V)\n",
    "    with vr_fun.x.petsc_vec.localForm() as vr_fun_local, \\\n",
    "            multiphenicsx.fem.petsc.VecSubVectorWrapper(vr, V.dofmap) as vr_wrapper:\n",
    "        vr_fun_local[:] = vr_wrapper\n",
    "    # Return eigenvalue and eigenfunction\n",
    "    return r, vr_fun"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78cdcdbb",
   "metadata": {},
   "source": [
    "In order to test this approach, recall that it can be shown by separation of variables in polar coordinates that\n",
    "$$\\eta^{(1)} = j_{0,1}^2$$\n",
    "where $j_{n, k}$ is the $k$-th positive zero of the $n$-th Bessel function $J_n$. We store the value of $j_{0,1}$ for later use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b674be5",
   "metadata": {},
   "outputs": [],
   "source": [
    "j01,  = scipy.special.jn_zeros(0, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f449b32",
   "metadata": {},
   "outputs": [],
   "source": [
    "j01, j01**2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aaa59794",
   "metadata": {},
   "source": [
    "Once $\\eta^{(1)}$ is known, the eigenvalue associated to any $\\alpha$ can be easily obtained as\n",
    "$$ \\eta^{(\\alpha)} = \\alpha j_{0,1}^2 $$\n",
    "by linearity."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a095ac72",
   "metadata": {},
   "source": [
    "We first test this approach by computing the eigenvalue for the case $\\alpha=0.1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ede319e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "eigenvalue_0p1, eigenfunction_0p1 = get_smallest_eigenvalue_dolfinx(0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5018fe3",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Computed: {eigenvalue_0p1}, vs expected {0.1 * j01**2}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80e5486a",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.isclose(eigenvalue_0p1, 0.578561)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88b60e9a",
   "metadata": {},
   "source": [
    "The approximation of $\\eta^{(0.1)}$ is reasonably accurate. We can also plot the corresponding eigenfunction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc38b676",
   "metadata": {},
   "outputs": [],
   "source": [
    "viskex.dolfinx.plot_scalar_field(eigenfunction_0p1, \"eigenfunction 0.1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9e0b8ff",
   "metadata": {},
   "source": [
    "We next try the case $\\alpha = 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4abed4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "eigenvalue_1, eigenfunction_1 = get_smallest_eigenvalue_dolfinx(1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d89dd67",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Computed: {eigenvalue_1}, vs expected {j01**2}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa6b49e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.isclose(eigenvalue_1, 1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9e40692",
   "metadata": {},
   "source": [
    "The current approximation of $\\eta^{(1)}$ is completely inaccurate. From the plot of the corresponding eigenfunction, we realize that the computed eigenvalue is a spurious one, which does not even correspond to an eigenfunction that satisfies the boundary conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "507b6e22",
   "metadata": {},
   "outputs": [],
   "source": [
    "viskex.dolfinx.plot_scalar_field(eigenfunction_1, \"eigenfunction 1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca49062c",
   "metadata": {},
   "outputs": [],
   "source": [
    "boundary_integral_1 = mesh.comm.allreduce(\n",
    "    dolfinx.fem.assemble_scalar(dolfinx.fem.form(eigenfunction_1 * ufl.ds)), op=mpi4py.MPI.SUM)\n",
    "print(f\"Boundary integral of the eigenfunction: {boundary_integral_1}\")\n",
    "assert not np.isclose(boundary_integral_1, 0.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3c233c1",
   "metadata": {},
   "source": [
    "For comparison, the boundary integral of the eigenfunction associated to $\\eta^{(0.1)}$ was indeed numerically zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c01255dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "boundary_integral_0p1 = mesh.comm.allreduce(\n",
    "    dolfinx.fem.assemble_scalar(dolfinx.fem.form(eigenfunction_0p1 * ufl.ds)), op=mpi4py.MPI.SUM)\n",
    "print(f\"Boundary integral of the eigenfunction: {boundary_integral_0p1}\")\n",
    "assert np.isclose(boundary_integral_0p1, 0.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2531d8d6",
   "metadata": {},
   "source": [
    "### Approach #1: fixup\n",
    "\n",
    "In the case $\\alpha = 1$, the previous computation shows a spurious eigenvalue. The value resulting from there is a consequence of the way Dirichlet boundary conditions are internally applied by `dolfinx` in `A = dolfinx.fem.petsc.assemble_matrix(..., bcs=[bc])`: the row/columns of `A` associated to every DOF on $\\partial\\Omega$ are cleared up, and the diagonal is set to 1. By repeating the same procedure the matrix `B` on the right-hand side, we end up introducing a spurious eigenvalue equal to 1.\n",
    "\n",
    "To obtain correct results, we should instead set the diagonal value of $A$ for row/columns associated to boundary DOFs to a number $d \\in \\mathbb{R}^+$, while still keeping a diagonal value equal to 1 for `B`. This will still introduce a spurious eigenvalue, with value $d$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e490c8dd",
   "metadata": {},
   "source": [
    "Let us try again for $\\alpha = 1$, setting $d = 1.5$. We then expect the computed eigenvalue to be equal to $d$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "354d3f79",
   "metadata": {},
   "outputs": [],
   "source": [
    "eigenvalue_1_fixup1, _ = get_smallest_eigenvalue_dolfinx(1.0, diagonal_A=1.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e046097d",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Computed: {eigenvalue_1_fixup1}, vs expected {j01**2}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcfdc061",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.isclose(eigenvalue_1_fixup1, 1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13f46fc4",
   "metadata": {},
   "source": [
    "We have successfully moved the first eigenvalue to $1.5$. Since we expect the true eigenvalue $\\eta^{(1)}$ to be around 5.78, the goal will be to choose $d$ large enough, so that the spurious eigenvalue does not get returned as the minimum one. For istance, choose $d = 10$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1defda0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "eigenvalue_1_fixup2, eigenfunction_1_fixup2 = get_smallest_eigenvalue_dolfinx(1.0, diagonal_A=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afd17896",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Computed: {eigenvalue_1_fixup2}, vs expected {j01**2}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e913e60",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.isclose(eigenvalue_1_fixup2, 5.785609)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69e8bcbc",
   "metadata": {},
   "source": [
    "A correct approximation of $\\eta^{(1)}$ is now obtained. Furthermore, due to linear algebra properties we expect the first eigenfunction of the case $\\alpha = 0.1$ and $\\alpha = 1$ to be the same, up to a sign."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67d877aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "viskex.dolfinx.plot_scalar_field(eigenfunction_1_fixup2, \"eigenfunction 1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c28404ce",
   "metadata": {},
   "source": [
    "In order to compare the eigenfunctions associated to the two different values of $\\alpha$, the next function normalizes eigenfunctions to ensure a consistent sign."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b59c7d2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def normalize(u: dolfinx.fem.Function) -> None:\n",
    "    \"\"\"Normalize an eigenvector.\"\"\"\n",
    "    scaling_operations: list[tuple[  # type: ignore[no-any-unimported]\n",
    "        dolfinx.fem.Function, typing.Callable[[dolfinx.fem.Function], ufl.Form],\n",
    "        typing.Callable[[petsc4py.PETSc.ScalarType], petsc4py.PETSc.ScalarType]\n",
    "    ]] = [\n",
    "        # Scale functions with a W^{1,1} norm to take away possible sign differences.\n",
    "        (u, lambda u: (u + u.dx(0) + u.dx(1)) * ufl.dx, lambda x: x),\n",
    "        # Normalize functions with a H^1 norm.\n",
    "        (u, lambda u: (ufl.inner(u, u) + ufl.inner(ufl.grad(u), ufl.grad(u))) * ufl.dx, lambda x: np.sqrt(x))\n",
    "    ]\n",
    "    for (function, bilinear_form, postprocess) in scaling_operations:\n",
    "        scalar = postprocess(mesh.comm.allreduce(\n",
    "            dolfinx.fem.assemble_scalar(dolfinx.fem.form(bilinear_form(function))), op=mpi4py.MPI.SUM))\n",
    "        function.x.petsc_vec.scale(1. / scalar)\n",
    "        function.x.petsc_vec.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7a49867",
   "metadata": {},
   "outputs": [],
   "source": [
    "normalize(eigenfunction_0p1)\n",
    "normalize(eigenfunction_1_fixup2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af34ae30",
   "metadata": {},
   "outputs": [],
   "source": [
    "viskex.dolfinx.plot_scalar_field(eigenfunction_0p1, \"eigenfunction 0.1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46030c6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "viskex.dolfinx.plot_scalar_field(eigenfunction_1_fixup2, \"eigenfunction 1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bf97385",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_errors(u: dolfinx.fem.Function, uu: dolfinx.fem.Function) -> None:\n",
    "    \"\"\"Compute errors between two different cases.\"\"\"\n",
    "    u_norm_form = (ufl.inner(u, u) + ufl.inner(ufl.grad(u), ufl.grad(u))) * ufl.dx\n",
    "    u_norm = np.sqrt(mesh.comm.allreduce(\n",
    "        dolfinx.fem.assemble_scalar(dolfinx.fem.form(u_norm_form)), op=mpi4py.MPI.SUM))\n",
    "    err_norm_form = (ufl.inner(u - uu, u - uu) + ufl.inner(ufl.grad(u - uu), ufl.grad(u - uu))) * ufl.dx\n",
    "    err_norm = np.sqrt(mesh.comm.allreduce(\n",
    "        dolfinx.fem.assemble_scalar(dolfinx.fem.form(err_norm_form)), op=mpi4py.MPI.SUM))\n",
    "    print(\"Relative error is equal to\", err_norm / u_norm)\n",
    "    assert np.isclose(err_norm / u_norm, 0., atol=1.e-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbcbd827",
   "metadata": {},
   "outputs": [],
   "source": [
    "compute_errors(eigenfunction_0p1, eigenfunction_1_fixup2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3464c3e2",
   "metadata": {},
   "source": [
    "### Approach #2\n",
    "\n",
    "While the fix introduced above is attractive because of its simplicity, it has a potential drawback: in cases where no analytical solution exists, the user has to try different values of $d$ until they found a value which is suitable for shifting away the spurious eigenvalues. \n",
    "\n",
    "In a second approach we leverage `multiphenicsx` capabilities to throw away DOFs associated to Dirichlet boundary conditions while assemblying `A` and `B`, so that they do not interfere with the eigenvalue calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35c97790",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Auxiliary function for the computation of the first eigenvalue using multiphenicsx\n",
    "def get_smallest_eigenvalue_multiphenicsx(  # type: ignore[no-any-unimported]\n",
    "    alpha: petsc4py.PETSc.RealType\n",
    ") -> tuple[\n",
    "    petsc4py.PETSc.RealType, dolfinx.fem.Function\n",
    "]:\n",
    "    \"\"\"Get the smallest eigenvalue, and the corresponding eigenfunction, using multiphenicsx.\"\"\"\n",
    "    # Define restrictions.\n",
    "    dofs_V = np.arange(0, V.dofmap.index_map.size_local + V.dofmap.index_map.num_ghosts)\n",
    "    dofs_V_partial_Omega = dolfinx.fem.locate_dofs_topological(V, boundaries.dim, facets_partial_Omega)\n",
    "    restriction_V = multiphenicsx.fem.DofMapRestriction(V.dofmap, np.setdiff1d(dofs_V, dofs_V_partial_Omega))\n",
    "    # Define problem forms\n",
    "    alpha_constant = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(alpha))\n",
    "    a = alpha_constant * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx\n",
    "    b = ufl.inner(u, v) * ufl.dx\n",
    "    # Assemble lhs and rhs matrices\n",
    "    A = multiphenicsx.fem.petsc.assemble_matrix(\n",
    "        dolfinx.fem.form(a), bcs=[], restriction=(restriction_V, restriction_V))\n",
    "    A.assemble()\n",
    "    B = multiphenicsx.fem.petsc.assemble_matrix(\n",
    "        dolfinx.fem.form(b), bcs=[], restriction=(restriction_V, restriction_V))\n",
    "    B.assemble()\n",
    "    # Solve\n",
    "    eps = slepc4py.SLEPc.EPS().create(mesh.comm)\n",
    "    eps.setOperators(A, B)\n",
    "    eps.setProblemType(slepc4py.SLEPc.EPS.ProblemType.GHEP)\n",
    "    eps.setDimensions(1, petsc4py.PETSc.DECIDE, petsc4py.PETSc.DECIDE)\n",
    "    eps.setWhichEigenpairs(slepc4py.SLEPc.EPS.Which.SMALLEST_REAL)\n",
    "    eps.getST().getKSP().setType(\"preonly\")\n",
    "    eps.getST().getKSP().getPC().setType(\"lu\")\n",
    "    eps.getST().getKSP().getPC().setFactorSolverType(\"mumps\")\n",
    "    eps.solve()\n",
    "    assert eps.getConverged() >= 1\n",
    "    # Extract first eigenvalue and eigenvector\n",
    "    vr = dolfinx.la.create_petsc_vector(restriction_V.index_map, restriction_V.index_map_bs)\n",
    "    vi = dolfinx.la.create_petsc_vector(restriction_V.index_map, restriction_V.index_map_bs)\n",
    "    eigv = eps.getEigenpair(0, vr, vi)\n",
    "    r, i = eigv.real, eigv.imag\n",
    "    assert abs(i) < 1.e-10\n",
    "    assert r > 0., \"r = \" + str(r) + \" is not positive\"\n",
    "    # Destroy EPS object\n",
    "    eps.destroy()\n",
    "    # Transform eigenvector into an eigenfunction that can be plotted\n",
    "    vr.ghostUpdate(addv=petsc4py.PETSc.InsertMode.INSERT, mode=petsc4py.PETSc.ScatterMode.FORWARD)\n",
    "    vr_fun = dolfinx.fem.Function(V)\n",
    "    with vr_fun.x.petsc_vec.localForm() as vr_fun_local, \\\n",
    "            multiphenicsx.fem.petsc.VecSubVectorWrapper(vr, V.dofmap, restriction_V) as vr_wrapper:\n",
    "        vr_fun_local[:] = vr_wrapper\n",
    "    # Return eigenvalue and eigenfunction\n",
    "    return r, vr_fun"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "447e921c",
   "metadata": {},
   "source": [
    "We test this approach on both cases $\\alpha = 0.1$ and $\\alpha = 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9c413be",
   "metadata": {},
   "outputs": [],
   "source": [
    "eigenvalue_0p1_approach2, eigenfunction_0p1_approach2 = get_smallest_eigenvalue_multiphenicsx(0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd99b90f",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Computed: {eigenvalue_0p1}, vs expected {0.1 * j01**2}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c7622eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.isclose(eigenvalue_0p1, 0.578561)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58d71db4",
   "metadata": {},
   "outputs": [],
   "source": [
    "eigenvalue_1_approach2, eigenfunction_1_approach2 = get_smallest_eigenvalue_multiphenicsx(1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c196f542",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Computed: {eigenvalue_1_approach2}, vs expected {j01**2}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04ff50e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.isclose(eigenvalue_1_approach2, 5.785609)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2b84633",
   "metadata": {},
   "source": [
    "We see that the computed values match the expected ones for $\\eta^{(0.1)}$ and $\\eta^{(1)}$. We also compare eigenfunctions between the two approaches."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a37ba051",
   "metadata": {},
   "outputs": [],
   "source": [
    "normalize(eigenfunction_0p1_approach2)\n",
    "normalize(eigenfunction_1_approach2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82a314a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "viskex.dolfinx.plot_scalar_field(eigenfunction_0p1_approach2, \"eigenfunction 0.1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a3704e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "viskex.dolfinx.plot_scalar_field(eigenfunction_1_approach2, \"eigenfunction 1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc938f36",
   "metadata": {},
   "outputs": [],
   "source": [
    "compute_errors(eigenfunction_0p1, eigenfunction_0p1_approach2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1f2581b",
   "metadata": {},
   "outputs": [],
   "source": [
    "compute_errors(eigenfunction_0p1, eigenfunction_1_approach2)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython"
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
